################### PROJECT- JCR REPLICATION - figures_and_table_main_01 #######################################
# This R file plots Figure 2, 3, 4, 5, 6, and Table 1 in the main article
# Those figures and table serve the descriptive purpose
# Below are brief introductions for those figures and table 
# Figure 2 - The spatial distribution of shifts in territorial control
# Figure 3 - The number of barangays with a shift in territorial control by year, between 2012 and 2014
# Figure 4 - The spatial distribution of battle-related violence
# Figure 5 - The number of barangays that experienced battle-related violence by year, between 2012 and 2014 
# Figure 6 - The spatial distribution of natural disasters
# Table 1 - The Distribution of Territorial Settings
# Last updated: 16th May, 2025, by Wangyin Zhao
##############################################################################################################

##############################################################################################################
# Outline of Script
# A. Preparation
#    - Set working directory
#    - Load required packages
# B. Import Data
# C. Figure 2
# D. Figure 3
# E. Figure 4
# F. Figure 5
# G. Figure 6
# H. Table 1
##############################################################################################################


#### A. Preparation ####
# Set working directory to Replication_Files folder
setwd("./Replication_Files")

# Import relevant packages
packages <- c("tidyverse", "sf", "gridExtra", "grid", "ggthemes", "knitr")
lapply(packages, library, character.only = TRUE)
########

#### B. Import data ####
# To generate those figures and table, the first step is to run the data processing files and to generate the analysis data 
# Upload dataset1_oldschool.csv for figure 2, 3, and 6 and Table 1
dataset1 <-  read_csv("./data/analysis/dataset1_oldschool.csv")

# Upload dataset2_oldschool.csv for figure 4 and 5
dataset2 <-  read_csv("./data/analysis/dataset2_oldschool.csv")

# The spatial data on the barangays from the Philippine Standard Geographic Code (PSGC) 
barangay <- st_read("./data/raw/Barangays/Barangays.shp") %>% st_make_valid()

# # The spatial data on the region from the Philippine Standard Geographic Code (PSGC) 
region <- st_read("./data/raw/regions/gadm41_PHL_1.shp")

########

#### C. Figure 2 ####
# 1. Data processing before plotting 
figure2_data <- dataset1 %>% 
  dplyr::select(ADM4_PCODE, year, Tchange_inf2) %>%
  filter(year != "2011") %>% 
  group_by(ADM4_PCODE) %>% summarize(sum = sum(Tchange_inf2))

figure2_data <- figure2_data %>%
  left_join(barangay) %>%
  # add the spatial geometry back 
  mutate(sum = as.character(ifelse(sum > 0, 1, 0))) %>%
  st_as_sf()

# 2. Plot figure 2
figure2 <- ggplot(figure2_data) +
  geom_sf(aes(fill = sum), size = 0.2, color = "grey", alpha = 1) +
  scale_fill_manual(values = c("grey", "blue"), name = "Shift in Territorial Control", labels = c("No shift", "Shift occurred")) +
  theme_minimal() +
  coord_sf(xlim = c(116, 127), ylim = c(4.3, 21)) +
  theme(legend.position = c(0.05, 0.94),  # Adjust the x and y coordinates
    legend.justification = c(0, 1),
    legend.title = element_text(size = 10))
########

#### D. Figure 3 ####
# 1. Prepare for the data 
figure3_data <- dataset1 %>% 
  dplyr::select(ADM4_PCODE, year, Tchange_inf2) %>%
  filter(year != "2011") %>% 
  group_by(year) %>% summarize(sum = sum(Tchange_inf2))

# Plot the figure
figure3 <- ggplot(figure3_data, aes(x = factor(year), y = sum)) + 
  theme_tufte(base_size=14) +
  geom_bar(stat = "identity", width = 0.25, fill = "grey50") +
  theme(axis.title = element_blank()) +
  geom_hline(yintercept = seq(0,800,200), col="white", lwd=1) +
  coord_flip() +
  theme(aspect.ratio = 0.5) +
  annotate("text", y=c(840,200,400,850), x=c(2.7,1.7,0.7,1.5), adj=1, family="serif",
           label=c("851","186","401","The number of \nthe shift of territorial control"))#

########

#### E. Figure 4 ####

# 1 Prepare for the data for the Davao region 
map <- barangay %>% mutate(ADM4_PCODE = as.integer(str_replace(ADM4_PCODE, "^PH0*", ""))) %>% dplyr::select(ADM4_PCODE) 
figure4_data <- dataset2 %>% filter(year != 2011) 

figure4a_data <- figure4_data %>% dplyr::select(ADM4_PCODE, year, npa_count) %>% 
  group_by(ADM4_PCODE) %>% 
  summarize(sum = sum(npa_count)) %>% 
  left_join(map)

figure4a_data <- figure4a_data %>%
  mutate(sum = as.character(ifelse(sum > 0, 1, 0))) %>%
  st_as_sf()

# Plot the Dvao region 
figure4a <- ggplot(figure4a_data) +
  geom_sf(aes(fill = sum), size = 0.2, color = "lightgrey") +
  scale_fill_manual(values = c("grey", "blue"), name = "Legend", labels = c("No Conflict", "Conflict")) +
  theme_minimal() +
  coord_sf(xlim = c(125, 126.5), ylim = c(5.2, 8)) +
  theme(
    legend.position = c(0.7, 0.25),  # Adjust the x and y coordinates
    legend.justification = c(0, 1),
    legend.title = element_text(size = 10)) +
  scale_x_continuous(breaks = seq(125, 126.5, by = 0.5))  # Adjust the labels as needed)


# Plot the rest of the Philippibes 
regions_to_color <- c("Davao del Norte", "Davao del Sur", "Davao Oriental", "Compostela Valley")
highlighted_data <- region[region$NAME_1 %in% regions_to_color, ]


figure4b <- ggplot() +
  geom_sf(data = region, color = "grey") +
  coord_sf(xlim = c(116, 127), ylim = c(4.3, 21)) +
  geom_sf(data = highlighted_data, fill = "blue", color = "black") +
  geom_rect(data = data.frame(xmin = 125, xmax = 126.6, ymin = 5.5, ymax = 8.2),
            aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
            fill = NA, color = "black", size = 0.3) +
  geom_text(data = data.frame(x = 125.5, y = 5, label = "Davao Region"),
            aes(x = x, y = y, label = label),
            color = "black",
            size = 2,
            hjust = 0.5, vjust = 0.5) +
  labs(x = NULL, y = NULL) +
  theme_minimal()

figure4 <- grid.arrange(
  figure4a, figure4b,
  widths = c(1, 2))

grid.lines(x = c(0.315, 0.5), y = c(0.35, 0.75), gp = gpar(lty = 2, col = "black"))
grid.lines(x = c(0.315, 0.5), y = c(0.25, 0.25), gp = gpar(lty = 2, col = "black")) ## Add a line between the plots
#########

#### F. Figure 5 ####
# 1. Prepare the data 
figure5_data <- dataset2 %>% filter(year != 2011) %>%
  group_by(year) %>%
  summarise(sum = sum(npa_count))

# 2. Plot the figure 
ggplot(figure5_data, aes(x = factor(year), y = sum)) + 
  theme_tufte(base_size=14) +
  geom_bar(stat = "identity", width = 0.25, fill = "grey50") +
  theme(axis.title = element_blank()) +
  geom_hline(yintercept = seq(0,90,30), col="white", lwd=1) +
  coord_flip() +
  theme(aspect.ratio = 0.8) +
  annotate("text", y=c(110,110,40,110), x=c(2.7,1.7,0.7,1), adj=1, family="serif",
           label=c("109","107","38","The number of \nbattle-related violence"))#
########

#### G. Figure 6 ####
# Prepare the data 
figure6_data <- dataset1 %>% 
  dplyr::select(ADM4_PCODE, year, ND_1.6_bar) %>%
  group_by(ADM4_PCODE) %>% summarize(sum = sum(ND_1.6_bar, na.rm = TRUE))

figure6_data <- figure6_data %>% 
  left_join(barangay) %>%
  mutate(sum = factor(ifelse(sum > 0, 1, 0), levels = c(0, 1))) %>%
  st_as_sf()

#2 Plot the figure 
figure6 <- ggplot(figure6_data) +
  geom_sf(aes(fill = sum), size = 0.2, color = "grey", alpha = 1) +
  scale_fill_manual(values = c("grey", "blue"), name = "The Distribution of Natural Disasters", labels = c("No", "Yes")) +
  theme_minimal() +
  coord_sf(xlim = c(116, 127), ylim = c(4.3, 21)) +
  theme(
    legend.position = c(0.05, 0.94),  # Adjust the x and y coordinates
    legend.justification = c(0, 1),
    legend.title = element_text(size = 10))

########

#### H. Table 1 ####
# 1.Prepare the data
IV_TS <- dataset1 %>%
  dplyr::select(ADM4_PCODE, year, T_inf2) %>% group_by(ADM4_PCODE) %>%
  mutate(T_inf2 = dplyr::lag(T_inf2, n=1)) %>%
  drop_na(T_inf2)

prop_table <- prop.table(table(IV_TS$T_inf2, IV_TS$year), margin = 2) # Create a table of proportions

percent_prop_table <- prop_table * 100 # Convert the proportions to percentages (multiply by 100)

rounded_percent_prop_table <- round(percent_prop_table, digits = 2) # Round the percentages to two decimal places (adjust as needed)

row_sums <- rowSums(table(IV_TS$T_inf2, IV_TS$year)) # Calculate the row sums (total counts) for each year

col_sums <- colSums(table(IV_TS$T_inf2, IV_TS$year)) # Calculate the column sums (total counts) for each T_inf2 value

total_count <- sum(row_sums) # Calculate the total count for the entire table

row_proportions <- row_sums / total_count*100 # Calculate the row proportions (proportion of each year)

row_proportions <- round(row_proportions, 2) # Round the row proportions to four decimal places

prop_table_with_sums <- cbind(rounded_percent_prop_table, Row_Sum = row_proportions) # Add row and column sums as proportions to the table

colnames(prop_table_with_sums)[ncol(prop_table_with_sums)] <- "Total" # Rename the empty column for better clarity

count_table <- table(IV_TS$T_inf2, IV_TS$year)

combined_table <- cbind(
  Count_2012 = count_table[, "2012"],
  Proportion_2012 = paste(rounded_percent_prop_table[, "2012"], "%", sep = ""),
  Count_2013 = count_table[, "2013"],
  Proportion_2013 = paste(rounded_percent_prop_table[, "2013"], "%", sep = ""),
  Count_2014 = count_table[, "2014"],
  Proportion_2014 = paste(rounded_percent_prop_table[, "2014"], "%", sep = ""),
  Total_Count = row_sums,
  Total_Proportion = row_proportions) # Combine the percentage table with counts

kable(combined_table, format = "latex", align = "c", 
      caption = "The Distribution of Territorial Settings",
      col.names = c("Count 2012", "Proportion 2012", 
                    "Count 2013", "Proportion 2013", 
                    "Count 2014", "Proportion 2014",
                    "Total Count", "Total Proportion"), digits = 4) # Use kable to generate a markdown table with both proportions and counts
########
